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SIMIARY 

Previous attempts to predict the transient acoustic pressure pulse at long 
horizontal distances from large explosions in the atmosphere have adopted a 
model atmosphere bounded above by a halfspace of finite sound speed and have 
represented the waveform as a superposition of contributions from dispersiyely 
propagating guided modes. Certain modes at low frequencies decay ejqjonentially 
(leaking modes) with increasing propagation distance. The practice up to now 
has been to neglect the contributions from such modes in such frequency ranges. 
The lower frequeney cutoffs for such modes are extremely sensitive to the 
nature of the upper halfspace in contradiction to the reasonable supposition 
that energy ducted in the lower atmosphere should be insensitive to the 
assumed form of the upper halfspace. In the present paper the overall problem 
is reexamined with account taken of poles off the real axis and of branch line 
integrals in the general Integral governing the transient waveform. Perturba- 
tion techniques are described for the computation of the imaginary ordinate of 
the poles and numerical studies are described for a model atmosphere terminated 
by a halfspace with c = 478 m/sec above 125 km. For frequencies less than 
0.0125 rad/sec, the GR^ mode, for example, is found to have a frequency 
dependent amplitude decay of the order of 10“^ nepers/km. Examples of 
numerically synthesized transient waveforms are exhibited with and without the 
inclusion .of leaking modes. The inclusion of leaking modes results in wave- 
forms with a more marked beginning rather than a low-frequency oscillating 
precursor of gradually increasing amplitude. Also, the revised computations 
indicate that wave foims invariably begin with a pressure rise, a result 
supported by other theoretical considerations and by experimental data. 

INTRODUCTION 

One of the standard mathematical problems in acoustic wave propagation is 
that of predicting the acoustic field at large horizontal distances from a 
localized source in a medium whose properties vary with height only. This 
problem, as well as its counterpart in electromagnetic theory, has received 
considerable attention in the literature (ref. 1) , is reviewed extensively in 
various texts (refs. 2-7), and, for the most part, may be considered to be 
well understood. 

A typical formulation of the transient propagation problem (refs. '8,9) 
leads (at sufficiently large horizontal distance r) to an intermediate result 
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which expresses the acoustic pressure as a double Fourier integral over angular 
frequency to and horizontal wave number k so that 

p = S(r) Re {/ £Cw)e'^‘^^ / [Q/DCa3,k)]e^^’^dkd£o}. (1) 

O -00 

Here S(r) is a geometrical spreading factor, which is X/t/x for horizontally 
stratified media and 1/ [agSin(r/a^)] ^ if the earth's curvatijre (a^ = radius 
of earth) is approximately taken into account. The quantity f(w) is the 
Fourier transform of a time- dependent function that characterizes the source. 

Q is a function of receiver and source heights Zj. and Zg, respectively, as well 
as of w and k, and possibly of the horizontal direction of propagation if winds 
are included in the formulation. In any case, given z^ and Zg, Q should have 
no poles in the complex k-plane when w is real and positive. The denominator 
D{(o,k) (which is termed the eigenmode dispersion function) may be zero for 
certain values k^^Cw) of k. 

The k integration contour for Eq. (1) is chosen to lie along the real 
k-axis except where it skirts below or above poles which lie on the real axis 
(see Fig. la, where branch lines are identified by slash marks, poles are 
indicated by dots, and the k integration contour is marked by arrowheads that 
show the direction of integration). Let it Suffice here to say that the placing 
of branch cuts and the selection of the contour must be such that the expression 
for the acoustic pressure dies out at long distance as long as a small amount 
of damping is included in the formulation. The guided-mode description in the 
formulation arises when the contour is deformed [permissible because of Cauchy's 
theorem and of Jordan's lemma (ref. 10)] to one such as is sketched in Fig. lb. 
The poles indicated there above the initial contour are encircled in the 
coimterclockwise sense, and there are contour segments which encircle (also in 
the counterclockwise sense) each branch cut that lies above the real axis. 

The integrals around each pole are evaluated by Cauchy's residue theorem so 
that what remains is a sum of residue terms plus branch line integrals. Each 
residue term is considered to correspond to a particular guided mode of 
propagation. 

One approximation that was previously made in the guided-mode formulation 
was to neglect contributions from poles [i.e», the kj^(w)] which were located 
above the real k-axis (refs. 8,9). The thought behind this omission was that 
most of the contributions in the synthesis of waveforms for long propagation 
distances would come from poles which were on the real k-axis. Another 
approximation was that, for long distances, the contribution from branch line 
integrals could be neglected as well. Given these two approximations, the 
expression for the acoustic pressure in Eq. (1) can be approximated as follows: 

‘^Un 

p = E S(r) / cos[wt - k^(oj)r + (j)j^(a))] du, (2) 

^ “Ln 

where Ajj(w) and <}!j^((jJ) Sixe defined in terms of the magnitude and phase of the 
residues of the integrand in Eq. (1) and the kj^(a3) are the real roots for 
D(w,k) (which are numbered in some order with n = 1,2,3, etc.). It is under- 
stood that in Eq. (2), for any given n, kn(w) should be a continuous function 
of 03 between the limits o)lj^ (lower) and 03^^^ (upper) . With this understanding, 
it should be possible to evaluate the resultant integral over o) approximately 
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by the method of stationary phase or by some numerical method. 

In spite of the seeming plausibility of the above two approximations, 
there is a set of circumstances intrinsic to low-frequency infrasonic propa- 
gation for which they are not valid, even for distances of propagation of 
more than 10,000 km. It is these circumstances and their relation to the 
analytic synthesis of guided-mode atmospheric infrasonic waveforms that are 
of central interest in this discussion. 

SYMBOLS FREQUENCY USED 

defined in Eqs, (6) 
sound speed for upper halfspace 
sound speed as a function of height 
eigenmode dispersion function defined in Eq. (5) 
defined in Eq. (3) 
gravitational modes 

horizontal wave number and its imaginary and real parts, 
respectively 
ordered roots of D(o),k) 
acoustic pressure 

horizontal distance of propagation 

[1,1] and [1,2] elements of the transmission matrix [R] , 
respectively 
time 

phase velocity Co3,k) 

complex phase velocity obtained by first iteration with Eq. {8a) 
roots of R]^j^(w,v) and R ]^2 > respectively 

roots of D(w,v) 
height 

height of bottom of upper halfspace 

derivatives of Rqi and R 12 with respect to v, respectively, and 
evaluated at v^ and v^^, respectively 
ambient density 
angular frequency 

characteristic frequencies used in Eq. (3) 
cutoff point in the (w/v) -plane for a non-leaking mode. 

INFRASONIC MODES 

An atmospheric model that is frequently adopted in studies of infrasound 
is one in which the sound speed c{z) varies continuously with height z in some 
reasonably realistic manner up to a specified height zq- and is constant (value 
c-p) for all heights exceeding Zj ^see Fig. 2). Should winds be included in 
the formulation, the wind velocities are also assumed to be constant in the 
upper halfspace z > z^. It would seem reasonable to say that there is some 
choice in specifying the values for both Zj and c-p, even though the computa- 
tions of such factors as Q and D(a),k) in Eq. (1) become more , lengthy with 
increasing Zp. Whatever the choice of Zp, it would seem reasonable to choose 
cp to be c(zp) so that the sound-speed profile would then be continuous with 
height. Intuitively, it would also seem that if the source and receiver are 
both near the ground and if the energy actually reaching the receiver travels 
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via modes of propagation channeled primarily in the lower atmosphere, then the 
actual value of the integral in Eq. (1) would be somewhat insensitive to the 
choices of zj and c-p. This idea, however, remains to be justified in any 
rigorous sense. In typical calculations performed in the past, Z'p was taken 
as 225 km, and c-p was taken as the sound speed (=« 800 m/sec) at that altitude 
(ref. 8). 

The formulation leading to that version of Eq. (1) which is appropriate 
to infrasound for frequencies at which gravitational effects are important 
(corresponding to periods greater than one to five minutes) is based on the 
equations of fluid dynamics with the inclusion of gravitational body forces, 
the associated nearly exponential decrease of ambient density and pressure with 
height, and a localized energy source. When c^ is taken to be finite, the 
incorporation of gravitational effects in this formulation leads to a disper- 
sion relation for plane waves propagating in the upper halfspace which is 
(winds neglected) (refs. 8,9) 

^2 = -G = [w - w^]/c.p - [w - o)g]k /o) , (5) 

where the solution of the linearized equations of fluid dynamics for z > z-p 
is of the form 

xlc z 

p/ /p^ = (Constant) e e^^^ e ^ . (4) 

In these equations p is again the acoustic pressure, is ambient density, x 
is the horizontal space dimension, and k 2 is the vertical wave number (alter- 
natively written as iG for inhomogeneous plane waves) . and Wg are two 
characteristic frequencies (co^ > cop) for wave propagation in an isothermal 
atmosphere where ooyy = (Y/2)g/cY and o)g = (y - 1)^'^ g/c^ (g 9.8 m/sec^ is 
the acceleration due to gravity and y « 1.4 is the specific heat ratio for air) 
The values of k (positive and negative) at which G^ is zero turn out to be the 
branch points in the k integration in Eq, (1). The branch lines extend i 5 )wards 
and downwards from the positive and negative branch points, respectively 
(recall Fig. 1) . 

The eigenraode dispersion function D(w,k) in the case of atmospheric 
infrasound can be written in the general form (ref. 8) 

D(w,k) = Aj2Rjp - - Rp2^* 

In this expression, and Rj ^2 elements of a transmission matrix [R] . 

They depend on the atmospheric properties only in the altitude range zero to 
z™, and are independent of what is assumed for the upper halfspace. In general 
their determination requires numerical integration over height of two simul- 
taneous ordinary differential equations [termed the residual equations (refs. 
8,9,11)]. They do depend on oa and k (or, alternately, on o) and phase velocity 
V = 03/k) , but are free from branch cuts. The other parameters A ^2 a^id A,, 
depend on the properties of the upper halfspace, and on w and k. Ajj and Aj 2 
are given (winds excluded) as 

= gk^/o)^ - yg/[2c^] ; (6a) 
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Ai2 = 1 - Cjk^/iJ. (6b) 

It may be noted that, since every quantity in Eq, (5) (with the possible 
exception of G) is real when oi and k are real, the poles that lie on the real 
k-axis (recall that they are the real roots of D) must be in those regions of 
the (a),k) -plane [or, alternatively, the Cu),v) -plane] where G^ > 0. Since at 
heighjts above zj, the integrand of Eq. (1) divided by j/pg should vary with z 
as e T, there is no leakage of energy into the upper halfspace for those 
modes that correspond to the above poles. Such modes are termed fully ducted 
modes . Modes for which there is leakage of energy are termed leaking . If D 
is considered as a function of u) and phase velocity v, the locus of its real 
roots v(a)) (dispersion curves) has [as has been found by numerical computation 
with the program INFRASONIC WAVEFORMS (ref. 8)] the general form sketched in 
Fig. 3. The nomenclature for labeling the modes (GR for gravity. S' for sound) 
is due to Press and Harkrider (ref. 12). It may be noted from Eq. (3) that 
there are two "forbidden regions" (slashed in the figure) in the (w. v)-plane. 
Within these regions there are no real roots of the function D(o 3 ,v) because G 
is imaginary. The existence of the high-frequency upper "forbidden region" 
implies that the phase velocities for propagating modes are always less than 
the sound speed chosen for the upper half space. The low-frequency lower-phase- 
velocity "forbidden region" appears to be due to the incorporation of gravita- 
tional effects into the formulation. However, if Cj is allowed to approach 
infinity, the lower "forbidden region" disappears. Thus, it can be seen that 
the fully ducted GRq and GRj^ modes both have a low-frequency cutoff [wl in Eq. 
(2)] which depends on c^. In fact, th<=^ larger Cj becomes, the smaller this 
cutoff frequency becomes. 

At this point, there should appear to be the following paradoxes. Given 
that frequencies below o)g may be important for the synthesis of a waVeform, an 
apparently plausible computational scheme based on the reasoning leading to 
Eq. (2) will omit much of the information conveyed by such frequencies. Also, 
in spite of the plausible premise that energy ducted primarily in the lower 
atmosphere should be insensitive to the choice for C'p, it can be seen that 
this choice governs the cutoff frequencies for certain modes and that certain 
important frequency ranges could conceivably be omitted by a seemingly logical 
choice for Cp. The resolution of these paradoxes seems to lie in the nature 
of the approximations made in going from Eq. (1) to Eq. (2). The latter 
equation may not be as nearly correct as earlier presumed, and it may be 
necessary to include contributions from poles off the real axis as well as from 
the branch line integrals. Even for the case when the propagation distance r 
is very long, it may be that the imaginary parts of the complex horizontal wave 
numbers are so small that the magnitude of e^^^ in Eq. (1) is still not small 
compared to unity. In addition, a branch line integral may be appreciable in 
magnitude at large r if there is a pole relatively close to the associated 
branch cut. 

ROOTS OF THE DISPERSION FUNCTION 

In light of the paradoxes mentioned, it would be desirable to modify the 
solution represented by Eq. (2) so as to remove the apparent artificial low- 
frequency cutoffs of the GRq and GR^ modes. As a first step, the nature of 
the eigenmode dispersion function D in the vicinity of the dispersion curve for 
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a particular mode is examined. The curve of values of phase velocity v 

versus co for a given (n-th) mode is known for frequencies greater than the low 
cutoff frequency wl. Given this curve, analogous curves and v^(o)) can 

be found for values of the phase velocity oo/k at which the functions Rij^(a),v) 
and C^) , respectively, vanish. One characteristic of the 

curves Vjj(oi)), Vg(aj) , and V|^((ji)3 which has been checked numerically for o) > 

(see Fig. 4) is that, for a given mode of interest, these curves all lie 
substantially closer to one another than to the corresponding curves for a 
different mode. 

Given the definitions above of Vg^(w) and v^((jJ), the dispersion relation 
D = 0 for a single mode may be approximately expressed, through a simple expan- 
sion, as 

D ~ - V ■ t^ll G](3) (v - v^) = 0, (7) 

where a = dR^/dv, and 3 = dRl 2 /dv, evaluated at v = Vg and v^, respectively 
(for simplicity, D is considered here as a function of co and v = oj/k rather 
than of ( 1 ) and k) , The above equation may also be written in the form 

vCl) = v^ + - v^)X/[l-X], (8a) 

where 

X = (3/a)(A^^ + G)/A^2V 

Eq. (8a) may be considered as a starting point for an iterative solution which 

develops v in a power series in Vg - v^. With v - Vg as the zeroth iteration, 

the right hand side of Eq, (8a) can be evaluated for the value of v required 
for the next iteration, etc. This iterative procedure should converge provided 
that Va or Vj^ is not near a point at which G vanishes and provided that G in 
the vicinity of Vg or Vj^ is not such that the variable X is close to unity. 

Among Other limitations, the iterative scheme is inappropriate for those 
values of w in the immediate vicinity of w^. 

As an illustration of the perturbation technique, detailed plots (for the 
GRq and GR]^ modes) versus angular frequency are given in Fig. 5 of w/ko (top 

portion of the figure) which is the reciprocal of the real part of 1/vtl), and 

of kj (bottom portion) which is the imaginary part of (kj^ and kj are the 

real and imaginary parts of k, respectively), where v^^J is the result of first 
iteration for the phase velocity using Eqs. (8). Note that kj is zero above 
the corresponding cutoff frequencies. The values shown in Fig. 5 are appropri-. 
ate to the case of a U. S. Standard Atmosphere (ref. 8; see also Fig, 2) 
without winds which is terminated at a height of 125 km by an upper halfspace 
possessing a sound speed of 478 m/sec. For frequencies at which Vg is computed, 
the agreement between vC^) and Vj^ has proven to be excellent. The w/kj^ serve 
as approximate extensions of the dispersion curves down to frequencies near 
zero, thus enabling the computation of waveforms with leaking modes included. 

TRANSITION OF MODES FROM NON-LEAKING TO LEAKING 

A more precise approximation to D(o 3 ,v) in the vicinity of cutoff [i.e., 
near the point (wl,Vl)] reveals that a dispersion curve becomes tangential to 
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the line G = 0 at (co^/v^). For w < there is a very narrow gap in the 
frequency range in which there are no poles in the k- (or v-) plane correspond- 
ing to a given n-th mode. This gap is of the order 10"1^ rad/sec for the GRq 
mode and 10"^ rad/sec for the GRj^ mode. 

Since there is a gap in the range of frequeneies for which a pole 
(corresponding to a mode) may exist, it is evident that evaluation of the 
integral over k in Eq. (1) by merely including residues may be insufficient for 
certain frequencies . Thus it would seem appropriate to include a contribution 
from branch line integrals. However, there is a line of reasoning which 
demonstrates that all contributions from branch line integrals are insignifi- 
cant as previously assumed. Further details on this matter are provided in 
reference 13. 


EXAMPLE (HOUSATONIC) 

Values of w/kj^ and kj calculated by the perturbation techniques outlined 
above were used [with a revised version of INFRASONIC WAVEFORMS (ref. 14)] to 
compute waveforms for the case of signals observed at Berkeley, California, 
following the Housatonic detonation at Johnson Island on October 30, 1962. 

A comparison of theoretical and observed waveforms for this case is given by 
Pierce, Posey, and Iliff (ref. 9). This case also serves as the main example in 
the 1970 AFCRL report by Pierce and Posey (ref. 8), and is discussed by Posey 
(ref. 15) within the context of the theory of the Lamb edge mode. The model 
atmosphere assumed here (winds included) is the same as in Fig. 3-12 of refer- 
ence 8, except that in the present model the upper halfspace begins at 125 km 
rather than at 225 km. To avoid repeating tedious calculations of the kj for 
the GRq and GRj^ modes for this model atmosphere, it was assumed that the kj 
would be close iii value to those shown in Fig. 5. 

In Fig. 6, sets of plots for the Housatonic case are shown with and without 
leaking modes. The waveform that includes leaking modes is regarded as an 
improvement in that among other things, the spurious initial pressure drop shown 
in the original waveform is not present here. In Fig. 7 of reference 9 
observed and theoretical waveforms are shown for the Housatonic case. On the 
basis of the calculations described above, this figure was redrawn and is 

given here as Fig. 7, The only difference between the two figures lies in the 

central waveform. The false precursor is absent in the waveform shown in Fig. 

7, and the first peak to trough amplitude has been changed from 157 bar to 

170 bar (less than a 10% increase). The remainder of the central waveform is 

virtually unchanged. The discrepancy with the edge-mode synthesis has not been 
diminished and remains a topic for future study. 

CONCLUDING REMARKS 

It was shown in this paper that, for a model atmosphere in which the sound 
speed is constant above some arbitrarily large height, the GRg and GRj^ modes 
have low cutoff frequencies and are leaking below that height. Given these 
facts, perturbation techniques were provided for the computation of the imagi- 
nary and real parts kj and kj^, respeetively, of the horizontal wave numbers for 
these modes . Knowledge of the kj and kj^ then made it possible to include, in 
a synthesis of waveforms, contributions from the GRq and GRj^ modes at frequen- 
cies where these modes were leaking. Finally it was demonstrated that this 
inclusion yielded waveforms that were more realistic than before. 
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Figure 6.- Numerically synthesized waveforms (Housatonic) 
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Figure 7.- Observed and theoretical 
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